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Abstract 

We consider the use of multigrid methods to solve the Euler equations for 
hypersonic flow. We consider the steady state equations with a Runge-Kutta 
smoother based on the time accurate equations together with local time step- 
ping and residual smoothing. We examine the effect of the Runge-Kutta co- 
efficients on the convergence rate considering both damping characteristics 
and convection properties. We also show the importance of boundary con- 
ditions on the convergence rate for hypersonic flow. Also of importance are 
the switch between the second and fourth difference viscosity. Solutions are 
given for flow around a bump in a channel and flow around a biconic section. 
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Introduction. 

Central difference type explicit schemes are currently being applied on a regular basis 
in the solution of the Euler and Navier-Stokes equations see, e.g [17]. In these applications 
multigrid is used to accelerate the basic explicit schemes. However, in most cases the 
applications have been to flows in the transonic or low supersonic regime. Higher speeds 
have generally been treated without multigrid or using semi-coarsening [8, 10]. 

Multigrid methods were first developed for elliptic equations. These were later ex- 
tended to hyperbolic equations such as the fluid dynamic equations for subsonic and 
transonic flow [4, 7, 12]. For these cases the steady state retains many of the properties 
of an elliptic equation since the region of supersonic flow is limited. We shall show that 
with proper care the multigrid method still works for hypersonic flow. Gustafsson and 
Lotstedt [2] have pointed out that hyperbolic multigrid works by two different processes. 
For the long waves the advection process is most important and multigrid achieves its 
efficiency by allowing the use of larger time steps on coarser grids. Hence, it is important 
that the smoother use large time steps. However, for the shorter waves dissipation is 
more important and the efficiency of multigrid is based on principles similar to that for 
elliptic equations. In this study we consider a Runge-Kutta scheme [3] as the smoother 
for the multigrid method. The central differences are augmented by an artificial viscosity 
based on TVD principles [13]. A number of changes are made to the scheme to enable it 
to work for higher speed flows. We consider both unbounded domains and domains with 
solid surfaces at the boundaries. The main difficulties are not due to the basic multigrid 
method but rather to the imposition of boundary conditions. We first describe the basic 
Runge-Kutta method for the central difference scheme together with a description of the 
artificial viscosity. We then describe the theory of the multigrid acceleration for a simple 
advection equation. We finally present several examples to demonstrate our conclusions. 

Basic Scheme 

The basic elements of the scalar dissipation model considered in this paper were first 
introduced by Jameson, Schmidt, and Turkel [3] in conjunction with Runge-Kutta ex- 
plicit schemes. The space discretization is based on central differences with an additional 
artificial viscosity. This algorithm has been used by many investigators to numerically 
solve the Euler equations for a wide range of fluid dynamic applications. The same type 
of space discretization has been applied to alternating direction implicit (ADI) schemes 
[9] and LU factored implicit schemes [5]. In this section the basic scheme is briefly 
reviewed. 

Consider the Euler equations in the form 

( 1 ) Wit + fx + 9y = 0 , 

where W is the four-component vector of conserved variables, and f y g are the flux 
vectors. The independent variables are time t and Cartesian coordinates (x,t/). If (1) 
is transformed to arbitrary curvilinear coordinates ( = £(x,y) and rj = ry(x,y), then we 
obtain 

(2) (J~ 1 W) t + F t + G 7l = 0 
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where J 1 is the inverse transformation Jacobian, and 

F = fVr, ~ gx , „ G = gx£ - fy v 

In a cell-centered, finite- volume method, (1) is integrated over an elemental volume in 
the discretized computational domain, and J -1 is identified as the volume of the cell. 
Equation (2) can also be written as 

J-'Wt + AW ( + BWr, = 0 

where A and B are the flux Jacobian matrices defined by A = dF/dW and B = dG/dW. 

To advance the scheme in time we use a multistage scheme. A typical step of a 
Runge-Kutta approximation to (2.2) is 

(3) WW = W (0) ~oc k ^ + D Tt G( k - 1) - AD] , 

where and D v are spatial differencing operators, and AD represents the artificial 
dissipation terms. The derivatives of the fluxes are approximated by central differences. 
Hence, we need to evaluate F at (i+ j). Since the dependent variables are given at (», j) 
we need to average them to calculate the flux at the cell face. This averaging can be 
done in several ways. In the original code the conservative variables were averaged to 
find the variables at the cell face and then the flux was evaluated. Some care must be 
taken with the geometric terms so that a constant field yields a zero derivative. In fact 
for the two dimensional case 5 quantities were averaged, the four conservative variables 
and E + p . This was done because of enthalpy damping. Alternatively one can average 
density, momentum and total enthalpy. This might be more accurate since the total 
enthalpy is constant in the steady state. An alternative to this procedure is to average 
the fluxes rather than the dependent variables again with modifications for metric terms. 
Averaging the fluxes allows for a more accurate representation of the Rankine-Hugoniot 
jump condition in one dimension. In this study we average the fluxes but have not seen 
much difference between the various approaches. 

The dissipation terms are a blending of second and fourth differences. That is, 


(4) 

AD=(Dl + Dl-D}-D*)w, 

where 

D\W = V< .,) A <] 

(5) 

(6) 

D\W = V ( [(a< (+w £\ A ( V ( A ( ] W, h 


and are the standard forward and backward difference operators respectively 

associated with the ( direction. The variable scaling factor A is chosen as 

( 7 ) X U+±,j = \ [( A €)i,j + ( A f).+i,i] ’ 
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where is proportional to the spectral radius of the matrix A. The coefficients e^ 2 ) and 
are adapted to the flow and are defined as follows: 


(8) 

e S+\,i “ * (2) 

max(i/,_ liJ , v itj , Vi+ij, v i+2 j), 

(9) 


\p i+ i,i - 2pi,i + Pi-i, A 

V *,J — 

\ Pi+l,j “1" ^PiJ “1“ Pi-1, j 1 

(10) 

e (4) = 

max [o, (*<*> - £>, ,)] , 


where p is the pressure, and the quantities and are constants to be specified. 
The operators for the rj direction are defined in a similar manner. 

The second-difference dissipation term is nonlinear. Its purpose is to introduce an 
entropy-like condition and to suppress oscillations in the neighborhood of shocks. This 
term is small in the smooth portion of the flow field. The fourth-difference dissipation 
term is basically linear and is included to damp high-frequency modes and allow the 
scheme to approach a steady state. Only this term affects the linear stability of the 
scheme. Near shocks it is reduced to zero. 

For high speed flows the switch (9) is not very good and does not allow the multigrid 
to converge. Instead we consider a TVD variation of the switch [13] given by 


( 11 ) 


. = \Pj±i ~ 2 Pj + Pi-i I 
** bi+1 - Pj\ + I Pi ~ Pi- 1 1 + € ' 


= 1 / 2 . 


With this change and the factor 1/2 in front of the second-difference dissipation term, 
the scalar equation becomes first-order upwind near shocks. In the case of the original 
v we find that v ~ .05 near shock waves in transonic flows, e must be chosen carefully 
to prevent the switch from being activated by noise. 

We now no longer have a free parameter for the second-difference dissipation. We 
also replace (10) by 

(i2) = max [°. * (4) (l ' re< i\,i)] - 

We choose T = 2 so that the fourth-difference dissipation is cut off for v > 1/2. The only 
free parameter is the coefficient of the fourth-difference term. 

Several other changes were made to the code in addition to the change to a TVD 
switch. In the original algorithm the artificial viscosity for the energy equation was 
based on the total energy rather than the energy. For high speed flows we base the 
artificial viscosity on the energy so that in each equation the basic dependent variable 
is also used in the artificial viscosity. This is more in line with upwind codes. This has 
previously been used in central difference codes [1]. The algorithm no longer preserves 
a constant total enthalpy in the steady state but enthalpy damping is not useful for 
supersonic flows. In most cases the difference between the two approaches is small with 
each approach having its advantages. The original version seems to give slightly sharper 
shocks. 
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To facilitate the treatment of boundaries, two layers of ghosts cells are created near 
each boundary. These ghosts cells are calculated by some extrapolation technique and 
then central difference are used at all interior points for the fluxes and artificial viscosities. 
At the solid boundary the normal velocity is reflected antisymmetrically to all ghost cells 
while the tangential velocity and density is reflected symmetrically to all ghosts cells. 
The pressure is calculated by the normal momentum equation. Hence, no special logic 
is required for the artificial viscosity near the solid surface. This procedure was found to 
be superior to setting the using zeroth order extrapolation for the first differences and 
setting the third difference equal to zero in calculating the artificial viscosity at the wall. 
We have also found that the effect of the second difference (5) is much more important 
near the solid surface than the fourth difference (6). Hence, there is no harm at setting 
the third difference equal to zero near the wall. 

We finally discuss the choice of the Runge-Kutta parameters a* and the time step 
in (3). For the coarser meshes and a scalar equation the central difference scheme with 
the artificial viscosity reduces to a first order upwind scheme. For this scheme TVD 
parameters have been suggested by Shu [11]. Tai [15, 18] has suggested other parameters 
that are optimal for damping characteristics. We have generally followed the suggestions 
of Tai as detailed later though the differences between the two set of coefficients was 
small. On the fine mesh we have a difficulty since in the smooth regions the scheme 
is a central difference operator with a fourth difference viscosity. However, near shocks 
the switch reduces the scheme to a first order upwind scheme. Hence, we have chosen 
Runge-Kutta coefficients that are stable for both sets of schemes. For most of the cases 
described below we have used a three stage Runge-Kutta time stepping scheme. On the 
coarser meshes we first considered the coefficients c*x = 1/9, c*2 = 1/3 and c*3 = 1. This 
has the property that for CFL = 3 we have pure convection, i.e. the wave moves one 
cell each stage without any dissipation. We then reduce the CFL to 1.5 to introduce 
dissipation. This turns out to be the same parameters suggested by Shu in his TVD 
scheme. After some experimentation we now choose c*x = .148,0*2 = .4 and 0*3 = 1 with 
no residua] smoothing and CFL = 1.5 as suggested by Tai. On the finest mesh we choose 
ax = .2,0*2 = -55 and c* 3 = 1 with a constant residual smoothing parameter 0 = .2 and 
CFL = 2.5 . In all cases the artificial viscosity was updated at every stage. 

We have also experimented with a time step that depends on Uij as given in (11). In 
this case the time step is smaller near strong shocks. We have also increased the second 
order difference on the coarser meshes by introducing a dependence on v, see also [10]. 

Multigrid 

Although we are primarily interested in solving the Euler equations, the simple ad- 
vection equation provides useful clues to the role of multigrid for accelerating the con- 
vergence. We next describe the multigrid theory for both the one and two dimensional 
model problems. The one dimensional analysis is used to motivate the choice of Runge- 
Kutta parameters and time step. The two dimensional analysis indicates the central role 
of the boundary conditions. 
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One dimensional theory 

Each Runge Kutta iteration has the potential to both convect and damp the error. In 
a multigrid scheme, the damping of the higher frequency errors is an essential property of 
the smoother. Lower frequency errors are more efficiently damped on coarser grids. For 
a fixed time step the convergence of the multigrid method is improved as the damping 
properties of the smoother are improved. On the other hand, to get the most advantage 
from the convective properties of the Runge Kutta scheme, a large time step is needed. 
An increase in the time step is often at the expense of the efficiency of the smoothing. To 
quantify these observations, we consider the one dimensional scalar advection problem 

(13) u t + u x = 0. 

We first analyze the periodic boundary condition case, where the Runge Kutta multi- 
grid scheme is optimized by considering only the damping properties of the RK smoother 
- convection cannot accelerate convergence. After discretizing the spatial derivative us- 
ing a central difference on a uniform mesh and adding both first and third order artificial 
viscosities, see (5) and (6). We rewrite equation (3) as 

(14) U W = - akAtRl*- 1 ') 

where 

( R (*“ 1 )V - u »+i ~ **»-! r^A;r U, ' +1 ~ 2u » + Ui ~ 1 

( h 2Ax (Ax) 2 

. J 4 )/* \ 3 u <+2 “ 4u i- 1 + 6 Ui - 4u i+ i + Ui-i 

+e ^ (Al) (Zxy • 

The convergence properties of the time iteration (used to solve the periodic problem) 
can be studied by looking at the amplification factor. If u = e' 5x / Ax is a Fourier mode, 
then the symbol of the residual is 

(15) R(6) = [isinfl + 4e^ 2 ^sin 2 ^ + 16e( 4 )sin 4 0] 

Ax 

and the symbol for the complete, k stage Runge-Kutta iteration is 

(16) u n+1 = g(6)u n 
where 

g(0) = 1 + a k z + ackdk-iZ 2 + b a k a *_i • • • oqz* 

and 

z = -AtR(0). 

For low frequency errors, |0| << 1, the symbol of the residual is also small. In fact, 

z « -A (id + 4eW), A = 
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and therefore |<7(0)| is close to one. Thus, low frequency errors are reduced very slowly. 
For high frequency errors, the size of |5'(^)| depends on the Runge-Kutta parameters as 
well as on the Courant number. To use the RK time stepping scheme effectively as a 
multigrid smoother, the amplification factor should be small for all frequencies which 
cannot be represented on coarser levels, i.e. the high frequencies. Various strategies, 
based on this model problem, have been used to obtain parameters which result in good 
Runge-Kutta smoothers. The basic philosophy is to uniformly reduce the amplification 
factor over all high frequencies, 7t/2 < 0 < 7r. We note, however, that the choices of the 
alpha’s and the Courant number depend very strongly on the type and strength of the 
dissipation terms. Optimal choices for a three stage scheme for the first order upwind 
discretization (e^ 2 ^ = 1/2, eO) = 0), for instance, results in an unstable scheme when used 
for the third order discretization with e^ 2 ) = 0, e^ 4 ) = 1/16. 

For linear problems, the multigrid algorithm consists of both the RK smoother and a 
coarse grid correction. In a typical two grid scheme, the coarse grid correction involves 
the restriction of the residual to a coarse level, solving for the correction on the coarse 
level, and then interpolation of the correction back to the fine grid. Multiple level 
schemes replace the coarse grid solve by another RK step together with a coarse grid 
correction to the next coarser level. Using reasonably good RK parameters with suitable 
restriction, interpolation and coarse grid operators results in a one dimensional iterative 
method which converges uniformly well for all mesh sizes. Moreover, the convergence 
per multigrid cycle can be made arbitrarily good using a RK smoother with a sufficient 
number of stages. This is not necessarily the case for the two dimensional problem, as 
we explain below. 

Two dimensional theory 

We now consider the two dimensional advection problem, 

(17) u t + au x + bu y — 0, 

first with periodic boundary conditions. Although the one dimensional analysis of the 
RK scheme can be used to get stable 2-d RK schemes, the two dimensional coarse grid 
correction is fundamentally different than in one dimension. There are low frequency 
errors whose residuals are magnified on coarser grids. The easiest way to understand this 
is the following. Suppose the space derivatives have been discretized on a uniform grid 
with cells of size h x h as in the one dimensional case, a central difference plus artificial 
viscosity in the form of both second and fourth differences. Then the Fourier symbol, 
corresponding to a Fourier mode, e x6x / h e x< fo / h ( Q f the two dimensional residual is 

hRh(0, (f>) = i(asin0 + 6sin^) + 4e( 2 )(asin 2 0 + 6 sin 2 (f > ) + 16e( 4 )(asin 4 0 + b sin 4 <f>) 
which, for small 0 and <f>, is 

(18) » i(a0 + b<j > ) + 4 eW(a0 2 + bcf> 2 ) + 16 e< 4 )(a0 4 + b<j> 4 ). 
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On the coarser grid, with cells of size 2 h x 2h, the symbol of the residual for the same 
frequency is 


hR 2 h(0, 4>) = 


i(a sin 2 9 + b sin 2 <f>) + /^^h) ( 2 ) ^ S * n2 ^ ^ sin 2 

2 (2/i) 2 

+16A(2 + 


For small 9 and <j> this is approximately 

(19) « i(a6 + b<j>) + 8 eW(a8 2 + b<f> 2 ) + 128e( 4 >(a0 4 + b<f> 4 ). 


Several observations can be made by comparing these two symbols, one representing 
the residual on the fine grid and the other representing the residual of the same error 
on a coarse grid. The Runge Kutta smoother cannot efficiently damp the low frequency 
errors. At the same time, the restriction and the interpolation operators transfer the 
low frequencies quite accurately. Thus, a two level multigrid cycle, even assuming that 
the coarse grid equations are solved exactly, can only solve for the lowest frequencies 
to the extent to which they can be represented on the coarse grid. Comparing the two 
expressions (18) and (19), we see the following. If a9 -\- btf> > Ch for some positive 
constant, C, then the first term in both expressions is the dominant term, and the 
expressions are approximately equal. Thus the coarse grid ’sees’ the same residual, 
and can correct for it. However, if a9 + b<f> is zero, i.e., the error is constant along 
characteristics, the imaginary part of the symbol vanishes and the situation is quite 
different. If only first order viscosity (i.e., ^ 0, — 0) is used on both the fine 

and the coarse grids, then, for small 6 and <f > , the ratio of the size of the residual when 
computed on the fine grid to the size of the residual computed on the coarse grid is 

^ 4e( 2 )(a0 2 + b(j) 2 ) 1 

* 8e( 2 >(a0 2 + b P) = 2 

If only the third order viscosity is used on both fine and coarse grids (i.e., = 0, 

e( 4 ) ^ 0), then this ratio is only 

16e^ (a9 4 + b<f> 4 ) 1 

W 128e( 4 )(a# 4 + b<j> 4 ) ~ 8' 

Since it is common to use the higher order viscosity on fine grids, and the less expensive, 
lower order viscosity on coarse grids, we also consider the case where e^ 2 ^ = 0 and e^ 4 ) ^ 0 
on the fine grid and e ^ ^ 0 and e^ 4 ) = 0 on the coarse grid. In this case, the ratio is 

_ 16e( 4 )(afl 4 + ty 4 ) 

8e( 2 )(a<? 2 + b<f > 2 ) 

Thus the corresponding multigrid algorithms will only be able to reduce the errors which 
are constant along characteristics by the approximate amounts of 1/2, 7/8 and 1 (not 
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at all), depending on the type of viscosities used. In particular, we note that this means 
that, no matter how well the RK scheme is optimized, the convergence rate per multigrid 
cycle cannot be made arbitrarily small. 

For non-periodic boundary conditions, for example those used for inflow/outflow or at 
a solid body, the situation is not as bad. In this case, convection plays an important role 
in the convergence process. To see the effect of inflow/ outflow boundary conditions on 
the convergence rates of the multigrid, consider the two dimensional advection equation, 
(17) on the unit square, with a = b = 1 with inflow boundary conditions specified at 
the bottom and left sides. Table 1 shows the difference in the convergence rates of the 
periodic problem and the inflow/outflow problem for a two grid multigrid algorithm. 
The computations were done by explicitly forming all of the matrices involved and then 
finding the eigenvalue of greatest modulus. We used the same RK coefficients for all 
cases, with ai = .2, = .55, 0*3 = 1 and A = 1.5. The coarse grid equations are 

solved exactly. The spectral radii are given for one to three RK smoothing steps per 
multigrid cycle. Table la is with e^ 2 ) = 1/2 , e^ 4 ) = 0 on the fine grid and e^ 2 ) = 1/2 
, = 0 on the coarse grid, Table lb is with e^ 2 ^ = 0 , e^ 4 ) = 1/16 on the fine grid 

and e^ 2 ) = 0, e^ 4 ) = 1/16 on the coarse grid, and Table lc is with e^ 2 ^ = 0 , e^ 4 ) = 1/16 
on the fine grid and e^ 2 ) = 1/2 , e^ 4 ) = 0 on the coarse grid. We make the following 
observations. The convergence rates of the inflow/outflow problems are much better 
than for the corresponding periodic problems. Rather surprisingly, the numbers indicate 
that this is true even as the problem size is increased. This effect is much greater than 
can be attributed to just the convection on the fine grid, since this would just reduce the 
convergence rate by a factor of 1 — 0(h). Moreover, additional smoothing can be used 
to further reduce the convergence rate for the inflow/ outflow problem. 

Model problem vs. real problem 

To what extent does the model problem analysis predict the behavior for our real 
problem of interest? In the model problem analysis we have tried to model the multigrid 
cycle as closely as possible. For the Euler equations, we use a cell-centered FAS multigrid 
approach for non-linear equations to accelerate the convergence of the Runge-Kutta time 
stepping to the steady state solution. For linear problems such as our model advection 
problem, FAS is equivalent to the correction scheme. A cell-centered scheme was used 
in both cases, with equivalent projection and interpolation operators for the case of 
uniform grids. There are many differences, however, which are more difficult to model. 
We mention a few of them here. 

Since the choice of the RK parameters depends strongly on the amount of artificial 
dissipation, it is somewhat surprising that the model problem parameters work at all for 
the Euler equations. A scalar viscosity is used for the system, the same amount for each 
equation. Thus the effective sizes of and e^ 4 ) which would be appropriate in the model 
problem analysis actually lie in a range of values. Similarly, although we use local time 
stepping, the Courant number, A* is a scalar quantity determined by the whole system 
and used for all four equations. Thus, the model problem should also be analyzed for a 
whole range of A’s. The use of the nonlinear switching mechanism on the fine grid also 
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RK steps 

periodic 

inflow/ outflow 

8x8 

16 x 16 

32 x 32 

8x8 

16 x 16 

32 x 32 

1 

.567 

.591 

.607 

.371 

.378 

.381 

2 

.517 

.519 

.513 

.138 

.145 

.154 

3 

.498 

.515 

.511 

.068 

.089 

.096 


RK steps 

periodic 

inflow/ outflow 

8x8 

16 x 16 

32 x 32 

8x8 

16 x 16 

32 x 32 

1 

.888 

.919 

.933 

.685 

.704 

.709 

2 

.859 

.898 

.919 

.477 

.507 

.527 

3 

.838 

.890 

.900 

.333 

.364 

.399 


RK steps 

periodic 

inflow/outflow 

8x8 

16 x 16 

32 x 32 

8x8 

16 x 16 

32 x 32 

i 

.951 

.990 

.998 

.697 

.748 

.768 

2 

.929 

.988 

.997 

.496 

.528 

.547 

3 

.911 

.986 

.997 

.354 

.389 

.401 


Table 1: Comparisons of spectral radii for two grid multigrid, periodic vs. inflow/outflow 
a), first order discretization on both fine and coarse grid levels b). third order discretiza- 
tion on both fine and coarse grid levels c). third order discretization on fine and first 
order on coarse grid level 

complicates the correspondence between the artificial viscosities for the model and real 
problems. Thus, for smooth flows, the model problem using e = 0 on the fine grid and 
e( 4 ) = 0 on the coarse grids has similar behavior. The non-uniform mesh with varying 
aspect ratios of the cells and the variation in the alignment of the flow with the grid can 
be partially modeled by considering the 2-d advection equation (17) with varying ratio, 
a/b. 

Results 

In all cases we use a 3 stage formula as described above. A simple V cycle is used 
with smoothing on the way down to coarser meshes and simple interpolation on the 
way up to finer meshes. The standard formula is to do one smooth on the finest mesh 
and 2 smoothings are all coarser meshes, though other choices are described later. The 
coefficient of the third order viscosity is k = 1/16. 

We first consider an inflow-outflow case. We use a uniform Cartesian grid with 
a constant inflow at the left and lower boundaries and outflow at the right and top 
boundaries. The Mach number and angle of attack are used to specify the three incoming 
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i=4 max=.57E-02 



i= 4 max=. 19E+00 




i= 4 max=.94E-02 




Figure 1: Entropy in inflow/outflow case, Mach 10, 45 degree angle: a) first order 
viscosity on fine and all coarse levels b) third order viscosity on fine and all coarse levels 
c) switched viscosity on fine grid and first order viscosity on all coarse levels. 





b). 



c). V(l,2,2) d). V(l,10,10) 




Figure 2: Flow over a bump in a channel, Mach 10 
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characteristic variables at inflow while the outgoing variable is extrapolated from the 
interior. At outflow one characteristic is incoming and is specified while the other three 
are extrapolated. The steady state solution is a constant. The initial conditions are the 
steady state solution plus a ten percent random perturbation. In figure 1 we plot the 
entropy for a flow entering the domain at 45 degrees to the x axis and a speed of Mach 
10 after 4, 8 and 16 multigrid cycles. One can clearly see that the main wave travels 
along the characteristics which for the entropy is the stream line, with the convergence 
rates consistent with the theoretical results for the two dimensional advection equation 
with various artificial viscosities. The effect of the advection is clearly seen. 

In the second case we consider flow in a channel with a 4 percent bump on the lower 
wall. The grid is uniform in the x direction and stretched in the y direction. In figure 
2a we show the grid while in figure 2b we plot the isomach lines for a Mach 10 flow. 
In figure 2c we show the convergence rate for both a 64 x 32 grid and a 128 x 64 grid 
with the standard multigrid parameters. In figure 2d we plot the convergence rate when 
we do one iteration on the finest mesh and 10 iterations on the coarser meshes instead 
of the standard 2 iterations on the coarser mesh. We see that doing more iterations 
on the coarser meshes does not do any harm. On the contrary in figure 2d we have a 
convergence rate independent of mesh size which we did not have in figure 2c. 

In the third case we consider flow about a cone consisting of two sections. The mesh 
is 128 x 32 with a nonuniform grid. The inflow is at Mach 6 with 0 degree angle of 
attack. In figure 3a we plot the convergence history for this case. In figure 3b we plot 
the residual of the continuity equation when the standard multigrid parameters are used. 
We have plotted both the entire field and a closeup of the leading edge. As expected the 
largest contributions to the residual come from the bow shock and two weaker shocks that 
originate at the point where the cone has a discontinuous tangent. In figure 3c we perform 
only one smoothing on the coarser grids. There is now a noticeable residual ahead of the 
bow shock in the free stream region. Thus, this multigrid scheme has caused disturbances 
to move upstream while in figure 3b there are only minor perturbations ahead of the 
shock. Hence, the central difference scheme allows disturbances to propagate upstream 
but the net effect is negligible. Hence, if the central multigrid acceleration does not 
contain enough smoothing on all grids then it can increase the effect of this upstream 
perturbation. Hence, there may be an advantage to considering an upwind multigrid 
even for central difference schemes. 

Conclusions 

The efficiency of multigrid for the Euler equations is governed by the smoothing 
properties of the basic scheme perpendicular to the characteristics and to convection 
properties of the scheme along the characteristics. On the coarse grids we have used 
Runge-Kutta parameters that are optimal for first order upwind schemes. On the finest 
mesh we use parameters that are stable for both a first order upwind scheme and a central 
difference scheme with a fourth difference artificial viscosity. By using residual smoothing 
we can increase the time step while maintaining good smoothing properties. In general 
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Figure 3a 
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Figure 3: Flow over a cone, Mach 6 
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we have found that a small amount of residual smoothing increased the robustness of 
the code. However, further increasing the residual smoothing and increasing the time 
step did not seem to improve the convergence rate. For a central difference scheme with 
artificial viscosity the second difference viscosity is important for the shock treatment. 
The behavior of the multigrid scheme is more affected by the fourth difference viscosity. 
Hence, a TVD scheme by itself is not a good smoother for a multigrid scheme until extra 
terms are added that can damp high frequencies. 

We have found that for inflow-outflow boundaries that the multigrid central differ- 
ence scheme converges rapidly for all Mach numbers. We have indeed shown that the 
convergence rates for the semi-infinite case are better than for the periodic case. The 
difficulties and slow convergence rates for more difficult geometries are due to boundary 
treatment. Care must all be exercised in the treatment of the artificial viscosity near 
the boundary. For the inflow- outflow case computing many iterations of the smoother 
on the coarse grids increases the convergence rate but not by very much so that it is 
not efficient. One drawback of the multigrid method is that it can convect disturbances 
upstream even for hypersonic problems. This propagation against the stream can be 
more pronounced than that introduced by the central difference smoother. 
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